/* Simulated annealing applied to the traveling salesman problem
 * Code by Florian Speelman & Jannis Teunissen */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double **dm;	//Making the distance matrix a global variable is almost necessary because it is a 2D matrix

double xrand(void) {	//Returns random number between zero and one, should not equal one
	static double fac = 1. / (double)(RAND_MAX + 1.);
	return (double)rand() * fac;
}

int readInput(const char *fname){ //Reads the input file
	int nlines, i, j;
	FILE *infile = fopen(fname, "r");
	if(!infile){
		printf("An error occurred opening %s\n", fname);
		exit(1);
	}
	
	fscanf(infile, "%d\n", &nlines);
	dm = (double**)malloc(nlines * sizeof(double*));
	for(i = 0; i < nlines; i++) {
		dm[i] = (double*)malloc(nlines * sizeof(double));
	}
	
	for(i = 0; i < nlines; i++){
		for(j = 0; j < nlines; j++){
			if(!fscanf(infile, "%lf", &(dm[i][j]))){
				printf("Something is wrong around element (%d,%d) in %s\n", i+1, j+1, fname);
				exit(1);
			}
		}
	}
	return nlines;
}

struct annealParameters {	//Struct containing the parameters for the simulated annealing
	double startTemp;
	double stopTemp;
	double pmove;
	double prot;
	double coolingCoef;
	int maxTotal;
	int itsAdjustT;
	int coolingType;
};
typedef struct annealParameters apars;


struct intdoub {
		int a;
		double b;
	};

int cmpStruct(const void *x, const void *y) {	//Used for the randOrder function to sort an array
	struct intdoub *p,*q;
	p = (struct intdoub*)x;
	q = (struct intdoub*)y;
	if(p->b > q->b) {
		return 1;
	}
	if(p->b < q->b) {
		return -1;
	}
	return 0;
}

void randOrder(int *x, const int N) {	//Randomize the order of x[]
	int i;
	
	struct intdoub *rlist;
	rlist = (struct intdoub*) malloc(N * sizeof(struct intdoub));
	
	for(i = 0; i < N; i++){
		rlist[i].a = x[i];
		rlist[i].b = xrand();
	}
	
	qsort(rlist, N, sizeof(struct intdoub), cmpStruct);
	
	for(i = 0; i < N; i++) {
		x[i] = rlist[i].a;
	}
}

double calcLength(int *x, int N) {	//Calculate the length of the path from c[1] to c[N] and back
	int i;
	double length = 0;
	for(i = 0; i < N; i++) {
		length += dm[x[i]][x[i+1]];
	}
	return length;
}

int accept(double dl, double *T, int *total, apars *pars){	//Returns 1 if the new path is accepted
	static int ctr = 0;
	ctr++;
	if(ctr > pars->itsAdjustT) {	//Sometimes lower the temperature
		ctr = 0;
		switch(pars->coolingType) {
			case 0: 
			*T = pars->stopTemp + pars->startTemp * pow((1. - (double)(*total)/pars->maxTotal), pars->coolingCoef);
			break;
			case 1:
			*T *= pow(pars->stopTemp / pars->startTemp, (double)pars->itsAdjustT/(double)pars->maxTotal);
			break;
			default:
			printf("Invalid cooling schedule\n");
			exit(1);
			break;
		}
	}
	if(dl < 0){
		return 1;
	}
	if(exp(-dl / (*T)) > xrand()) {
		return 1;
	}
	return 0;
}

//Used to rotate a segment
void rotateSegment(int **cptr, int N, int start, int end){
	int n, y, *c;
	c = *cptr;
	for(n = 0; n < (end-start + 1)/2 ; n++) {
		y = c[start + n];
		c[start + n] = c[end - n];
		c[end - n] = y;
	}
	c[0] = c[N]; c[N+1] = c[1];	
}

//Used to move a segment to another place
void moveSegment(int **cptr, int **c2ptr, int N, int start, int end, int newPos){
	int n, m, *c, *c2;
	c = *cptr;
	c2 = *c2ptr;

	if(newPos < start){
		n = 1;
		for(m = 0; m <= end - start; m++){
			c2[n++] = c[start + m];
		}
		for(m = newPos + 1; m < start; m++){
			c2[n++] = c[m];
		}
		for(m = end + 1; m <= N; m++){
			c2[n++] = c[m];
		}
		for(m = 1; m <= newPos; m++){
			c2[n++] = c[m];
		}
	} else {
		for(n = 1; n < start; n++){
			c2[n] = c[n];
		}
		for(m = end + 1; m <= newPos; m++){
			c2[n++] = c[m];
		}
		for(m = 0; m <= end - start; m++){
			c2[n++] = c[start + m];
		}
		for(m = newPos + 1; m <= N; m++){
			c2[n++] = c[m];
		}
	}
	c2[0] = c2[N]; c2[N+1] = c2[1];	
	*cptr = c2;
	*c2ptr = c;
}



double anneal(apars *pars, int N, int **cityOrder) {	//The function that does the actual annealing
	int ctr, n, start, end, y, *c, *c2, newPos, *best;
	double r, length, minLength, T, dl; 
	
	T = pars->startTemp;
	
	c = (int*) malloc((N + 2) * sizeof(int));	//The sequence of cities is stored in c[]
	c2 = (int*) malloc((N + 2) * sizeof(int));	//Used as a copy of c[]
	best = (int*) malloc((N + 2) * sizeof(int));	//The best solution
	
	for(n = 1; n < N + 1; n++) {
		c[n] = n-1;	//Fill c[] with the number 1-N, c[0] and c[N] are to make accesing c[n+1]  and c[n-1] easier
	}
	
	randOrder(&c[1], N);	//Randomize the sequence in c[]
	c[0] = c[N]; c[N+1] = c[1];
	length = calcLength(c, N);
	minLength = length;
	ctr = 0;
	
	while(ctr < pars->maxTotal) {
		ctr++;
		
		do {	//generate a start and endpoint in the list
			start = (int)(xrand() * N + 1);
			end = (int)(xrand() * N + 1);
			if(start > end) {
				y = end;
				end = start;
				start = y;
			}
		} while (end - start > N * 0.5 || end == start);
		
		r = xrand(); //To choose the operation that generates a neighbouring state
		
		if((r -= pars->prot) < 0) { // Rotate section
			dl = dm[c[start-1]][c[end]] + dm[c[start]][c[end+1]] - dm[c[start-1]][c[start]] - dm[c[end]][c[end+1]];
			if(accept(dl, &T, &ctr, pars)) {
				length += dl;
				rotateSegment(&c, N, start, end);
			}
		} else if((r -= pars->pmove) < 0) { // Move section
			newPos = (int)(end + xrand() * (N - (end - start + 2))) % N + 1;
			dl = dm[c[newPos]][c[start]] + dm[c[end]][c[newPos+1]] + dm[c[start-1]][c[end+1]]
				- dm[c[start-1]][c[start]] - dm[c[end]][c[end+1]] - dm[c[newPos]][c[newPos+1]];
			if(accept(dl, &T, &ctr, pars)) {
				length += dl;
				moveSegment(&c, &c2, N, start, end, newPos);
			}
		} else { // Swap two cities in the list
			if(end - start == 1) {
				dl = dm[c[start-1]][c[end]] + dm[c[start]][c[end+1]] - dm[c[start-1]][c[start]] - dm[c[end]][c[end+1]];
			} else {
				dl = dm[c[end-1]][c[start]] + dm[c[start]][c[end+1]] + dm[c[start-1]][c[end]] + dm[c[end]][c[start+1]] 
					   - dm[c[start-1]][c[start]] - dm[c[start]][c[start+1]] - dm[c[end-1]][c[end]] - dm[c[end]][c[end+1]];
			}
			if(accept(dl, &T, &ctr, pars)) {
				length += dl;
				y = c[start];
				c[start] = c[end];
				c[end] = y;
				c[0] = c[N]; c[N+1] = c[1];
			}
		}
		
		if(length < minLength) {
			minLength = length;
			for(n = 0; n < N + 2; n++) {
				best[n] = c[n];
			}
		}		
	}
	free(c);
	free(c2);
	*cityOrder = best;
	return minLength;
}
		

int main(void) {
	int N, i;	//Number of places to visit
	int *cityOrder;	//The optimal ordering of the cities is stored in cityOrder
	double minLength;
	srand(1337);
	
	apars parameters;
	parameters.startTemp = 10000.;
	parameters.stopTemp = 0.1;
	parameters.maxTotal = 1000*1000;
	parameters.itsAdjustT = 10*1000;
	parameters.prot = 0.75;
	parameters.pmove = 0.25;
	parameters.coolingType = 1;
	parameters.coolingCoef = 4.;
	
	N = readInput("citydm.txt");
	minLength = anneal(&parameters, N, &cityOrder);	//Note that N can be set to a lower value manually for simulating less cities
	
	printf("The minimum distance found is %.2f km\n", minLength);
	printf("The ordering of the cities is as follows:\n");
	for(i = 1; i < N+1; i++){
		printf("%d -> ", cityOrder[i]);
		if(i % 10 == 0)puts("");
	}
	printf("and back to %d again...\n", cityOrder[1]);
	
	return 0;
}
		
	
